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Abstract 



^ ' We investigate the ability of the reference hypernetted-chain inte- 

B' 

gral equation to describe the phase diagram of square-weh fluids with 



four different ranges of attraction. Comparison of our results with 
simulation data shows that the theory is able to reproduce with fairly 
good accuracy a significant part of the coexistence curve, provided an 
I extrapolation procedure is used to circumvent the well-known patholo- 

CSI ' gies of the pseudo-spinodal line, which are more severe at reduced 

■ width of the attractive well. The method provides a useful approach 

. for a quick assessment of the location of the liquid-vapor coexistence 



curve in this kind of fiuid and serves as a check for the more complex 



H ' problem of anisotropic "patchy" square- well molecules. 
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Introduction 



Structural and thermodynamic properties of the square-well (SW) fluid have 
been studied with a huge variety of statistical mechanical methods. This 
system has long been viewed as the simplest nontrivial model able to capture 
the main phenomenology of real atomic fluids by complementing hard-sphere 
repulsion with a flnite-range and constant attractive well [T]. Almost all 
theories for the hquid state have been applied to the study of the SW fluid 
and many numerical simulation studies have been published [21 El HI El El [71 
[HI M [ini [HI [121 [131 [E] , transforming this system into an important testbed 
for theories [21 HSl HSl [13 [HI |19l 120] • 

More recently, the study of colloids and protein solutions has renewed 
interest in SW fluids as a tool to study trends and general features of strongly 
localized isotropic attractions. However, more detailed investigations of such 
systems suggest that the isotropic model should be modifled to account for 
highly directional (patchy) interactions [21]. The recent model studied by 
Kern and Frenkel [22] can be described as a SW model where the potential 
well has finite an^ii/ar extent as well as finite rac/m/ extent. 

A key problem with such patchy fiuids is the determination of their phase 
diagram. Specialized computer simulation methods have been developed to 
allow efficient and reliable determination of liquid- vapor coexistence. While 
applications to atomic systems is nowadays an almost trivial exercise, numer- 
ical simulation in the presence of highly attractive and directional forces is 
still a challenging and time-consuming task. Theoretical modeling, albeit ap- 
proximate, may provide a worthwhile tool for faster and wide-ranging scans 
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of the relevant parameter space. For such purposes, perturbation methods 
have been used extensively in the past, but their accuracy in the presence of 
strongly directional attractions is not uniform. Alternatively, modern inte- 
gral equation theories are known to provide reliable and accurate description 
of structural properties for isotropic [l7l [23] as well as anisotropic poten- 
tials [21] . It is true that in the vicinity of a liquid-vapor critical point they 
manifest shortcomings ranging from branching of multiple unphysical solu- 
tions [25j to wrong critical behavior [T7]. Notwithstanding such limitations, 
we believe that integral equation approach is still able to provide a quick 
and reasonably rehable description of the phase diagram. Although results 
for fluid-fluid coexistence from modern integral equation theories are rela- 
tively scarce, the existing evidence shows that they may provide an accurate 
description of the low temperature part of the binodal curve and, by ex- 
trapolation, make it possible to approximately locate the critical point in 
one-component p!7l [26] and two-component systems [27]. More accurate ap- 
proximations, such as SCOZA or HRT, have been applied to the SW fluid 
[T9l [20] , but their implementation for anisotropic interactions is not an easy 
task. 

Motivated by an interest in applying the reference hypernetted chain 
(RHNC) approximation |23J to simple anisotropic models of patchy colloids, 
as a preliminary step we wanted to understand the limits and shortcomings 
of the integral equation approach to such systems. For this reason, we have 
undertaken here the study of the liquid- vapor coexistence of an isotropic SW 
fluid, the extreme limit of patchy SW potentials, using the RHNC approxi- 
mation. 



2 



The RHNC approximation has been used before to study SW fluids [I6l 
but, so far as we know, it has not been used to study phase coexistence in 
such systems nor it has been tested at temperatures below the critical tem- 
perature. In this communication, we report results for structural properties 
and liquid- vapor coexistence of a few SW fluids of varying attraction range. 
The plan of the paper is the following: In section 2, we summarize the basic 
RHNC theory and discuss the computational details of our calculation of the 
phase coexistence curve. Results are presented and critically discussed in 
section 3. A short summary of our findings and conclusions is assembled in 
the final section. 



1 RHNC theory and phase coexistence cal- 
culations 

We consider a system of spherical particles of diameter a interacting via a 
pair-wise square-well potential 0(r) given by 



oo, r < a 



[r) = <-e, cr < r < Act, 



0, Xa <r 



where e is the depth and A the dimensionless extent of the potential well. 
Reduced temperature and density are introduced as usual using e and a as 
energy and length scale respectively: T* = kBT/e, p* = pa^. 
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1.1 RHNC integral equation 

The pair distribution function g{r) of a classical fluid is related to the pair 
potential 0(r) by the exact relation [28] 

g{r) = exp [-p(j){r) + h{r) - c(r) + B(r)] , (1) 

where /3 = l/ksT is the inverse Kelvin temperature, h{r) = g{r) — 1 is the 
pair correlation function, and c(r) is the direct correlation function defined 
via the Ornstein-Zernike (OZ) equation, 

h(r) = c(r) P J dr'c(|r — r'\)h{r'). (2) 

It is more convenient for numerical work to solve equations ([1]) and ([2]) for 
the indirect correlation function 7(r) = h{r) — c(r) with the OZ equation 
de-convoluted in Fourier space; the pair of equations to solve then becomes 

c(r) = exp [-p^ir) + 7(r) + B{r)] - 1 - 7(r), (3) 
These equations are connected through the transforms 
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c(^) = ^ / dr rc{r) sm{kr) (5) 



k Jo 
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following equation and 



7(r) 



1 




■oo 



dk k'^{k) sin(fcr) 



(6) 



27r2' 



following equation (jl]) to form an iteration cycle. The so-called bridge func- 
tion B{r) that appears above is a complicated functional of the pair cor- 
relation function for which no exact computable expression is known. The 
various integral equation formulations found in the literature correspond to 
different approximate forms for B{r). 

The SW direct correlation function c(r) is of course discontinuous at r = a 
and r = Act. To compute its Fourier transform ([5]), we assign the single value 
of c(r) at a discontinuity to be the arithmetic mean of its separate values at 
the discontinuity |2H|; e.g., 




We further ensure that the discontinuities fall on calculated grid points in 
r. Moreover, note that it is only the short-ranged direct correlation function 
c(r) that is expected to vanish at large r; i.e., its Fourier transform requires 
c(rmax) = 0, where rmax is the cutoff distance in the r grid. In the discrete 
notation used below, rmax = NrAr. The pair correlation function h{r), on 
the other hand, may have whatever long-range tail it wishes without affecting 
the calculation. 

The RHNC closure [23] assumes that the unknown bridge function B{r) 
can be approximated by the corresponding known bridge function BQ{r) of a 




c{Xa) = lim 

£— >0 



c{Xa + e) + c{Xa — e) 



(7) 
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reference system. Within this approach, just as in the original hypernetted- 
chain closure [30], the excess Helmholtz free energy per particle can be written 
in closed form as 



N N N N 



where 



pAi _ I f f 1 
f3A2 1 f dk 



N 2pJ (27r) 

PAs PAl 



-p J dr |-/^2(r) + h{r) - g{r) In [^7(r)e^^M] } , (9) 
J-^{\n[l + ph{k)]-ph{k)}, (10) 



N N 



-^pj dr[g{r)-go{r)]Bo{r). (11) 



Equation (ITT]) directly expresses the RHNC approximation; here A^ is the 
corresponding reference system contribution, computed from the known free 
energy A^ of the reference system as A^ = A^ — A'l — A2, with A^ and A2 
calculated as above but with reference system quantities. 

Length and energy parameters of the reference system, respectively ctq 
and eo, should be optimized by variational minimization to obtain the free 
energy from equation ([8]), which leads [H] to the conditions (written here in 
dimensionless forms) 

p I dv[g{r) - go{r)]ao^^ = 0, (12) 
p J dv[gir) - go{r)]eo^^ = 0. (13) 

These guarantee thermodynamic consistency between the direct calculations 
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of the pressure P and internal energy U using g{r) 



(5P = p 




I 



dr g{r)r 



(14) 



dr 



U_ 

N 




dr g{r)(f){r) 



(15) 



and the corresponding thermodynamic derivatives of the free energy, (3P — 
p = p^d{(3Ae^/N)/dp and U/N = d{f3Ac^/N)/df3. At present, however, only 
the hard sphere (HS) system, with sphere diameter ctq but no energy scale, 
is well enough known to serve as reference system. With this choice, only 
equation (|T^ is implemented here for the results of the next section. 

We will thus obtain the two main ingredients needed, pressure and chem- 
ical potential, in a consistent way with no additional approximations be- 
yond the choice of hard spheres as reference system. Any inaccuracy in the 
computed results can be ascribed solely to the quality of the chosen bridge 
function. 

We solve the RHNC equations numerically on r and k grids of A^^ = 2048 
points with intervals Ar = O.Olo" and Ak = n /{NrAr) using the combination 
of Newton-Raphson and Picard methods introduced by Gillan [32] and op- 
timized by Labfk et al. [33]. Fourier transforms are evaluated with the Fast 
Fourier Transform algorithm. Sample calculations with the same Ar = O.Olcr 
but increased A^^ = 4096 and A^ = 8192 produce identical results for all the 
printed output; i.e., to four or five significant figures. 

For the bridge function, we use the Verlet-Weis-Henderson-Grundke 
( VWHG) parametrization of HS numerical simulation results [Ml [35] . The 
HS free energy is taken from the Carnahan-Starling equation of state [36] , 
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which is also incorporated into the VWHG parametrization. Other parametriza- 
tions exist and in the only previous investigation of the SW fluid with RHNC, 
at supercritical temperatures, Gil-Villegas et al. used that of Malijevsky 
and Labi'k [371 EHl and found excellent agreement with computer simulation 
results for thermodynamics and correlation functions. In the present study 
we choose VWHG bridge functions because they have been frequently used 
in RHNC calculations (also for liquid-vapor coexistence of HS Yukawa [39] 
fluids) , thus allowing a direct comparison of the performances of the approx- 
imation with different potentials. 

1.2 Liquid-vapor coexistence 

Two-phase coexistence at constant temperature requires the equality of pres- 
sure and chemical potential of the two phases; RHNC provides convenient 
expressions for both quantities, 



where we have now specialized the pressure equation (HM for the SW poten- 
tial and we are using the cavity function y{r) = g{r) exp[/?0(r)], which has 
no discontinuities. An additional density-independent term from the ideal 
gas limit has been neglected in the total chemical potential fi. 

In order to find the densities of the two coexisting phases, we use the 
following procedure [41J: we start two calculations, one at low density and 




A3y(Ao-)(e^^-l) , 



(16) 




(17) 
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the other at high density at the temperature of interest. Then, density is 
increased from the lowest value and decreased from the highest until either 
the coexistence conditions are satisfied or the numerical program is unable 
to converge, even after systematic reduction of the density change, signaling 
the disappearance of a solution of the integral equation. Conditions of equal 
pressure and chemical potential are verified by looking for numerical solutions 
of the nonlinear set of equations, 

/3P(p,) = /?P(pg), (18) 
= Pfi{p,), (19) 

where Pg and pi denote respectively the density of the gas and of the liq- 
uid. The possibility of coexistence can be visually checked by plotting in a 
{PP, Pp) plane the high and the low density branches. A crossing of the two 
curves signals the occurrence of phase coexistence. The exact determina- 
tion of the intersection point is performed numerically by using cubic spline 
interpolations through the calculated points. 

Ideally, one should get crossings along the whole binodal line, from lowest 
accessible temperature up to the critical point. In practice, the approximate 
nature of RHNC, shared with all integral equations, introduces thermody- 
namic inconsistencies and modifies the singular behavior of the spinodal line 
obtained from the fluctuation route. It is transformed into a pseudo-spinodal 
line of branching points where the real solution of the equations does not show 
a divergence of the structure factor S{0) = 1 -|- ph{0) at zero wavenumber 
or, equivalently, a vanishing of the inverse long-wavelength (Iw) isothermal 
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compressibility 



P^^^^l-pc{0). (20) 



Moreover, such a pseudo-spinodal line extends the no-solution region into 
a non-negligible neighborhood of the critical point, thus preempting the pos- 
sibihty of finding the coexisting densities. As a result, the possibility of get- 
ting the coexistence curve with integral equations is limited to relatively low 
temperatures. In a way, the quality of an integral equation can be accurately 
tested by the extent of the resulting coexistence curve. 

In the case of SW fluids we have found, as discussed in the next section, 
that reducing the width of the SW progressively extends the temperature 
interval below the critical point where the pseudo-spinodal line preempts 
the binodal. However, we note that the quality of the thermodynamic and 
structural data in the accessible region of the phase diagram remains high. 
Thus, we have investigated the possibility of a small extrapolation of the 
RHNC data inside the pseudo-spinodal region. In all the cases we have 
investigated, such extrapolation is required only in the gas phase and the 
resulting coexistence curve appears as a smooth continuation of the part 
based on real crossings. 

In figure 1, we show a typical case of quasi-crossing while figure 2 shows 
a case of a more distant missed crossing. In both cases, the possibility of 
finding a coexisting vapor phase is preempted by the sudden disappearance 
of the low-density solution. It is clear that, although RHNC does not have a 
low-density solution at the same temperature and pressure of a corresponding 
liquid, a smooth extrapolation of thermodynamic data could be quite safe in 
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• gas 
■ liquid 



-4 1 ' ' ' ' ' ' ' ' ' ' ' ^ 

0.2 0.4 

Figure 1: Determination of the coexistence point. Case of quasi-crossing at 
T* = 2.5, A = 2. Squares: liquid branch; circles: vapor branch. Dashed 
lines indicate the cubic spline approximation used to locate the extrapolated 
crossing point. 



11 



1 1 1 1 1 1 1 1 1 1 



7 ' / 






/ 








• g^^ 
■ liquid 




1 , , , , 1 







I \ \ I \ I \ \ \ \ I \ \ \ \ I 
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Figure 2: Determination of the coexistence point. Case of a large extrapola- 
tion at T* = 1.74, A = 1.75. Symbols as in figure 1. 
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cases like that in figure 1, due to the small curvature of the lines and to their 
strong transversality. Larger extrapolations, like that in figure 4, appear to 
introduce uncontrolled uncertainties. Additional comments will be added in 
the next section. 

2 Results 

As a first check of the RHNC quality, we have computed the pressure P and 
excess chemical potential 



for states at high densities with A = 1.5 examined by Labik et al. [40j using 
scaled-particle Monte Carlo (SP-MC) simulation at p* = 0.8 and 0.9. Results 
for the pressure from simulation and the RHNC calculation are compared in 
figure 3, while figure 4 gives the comparisons for the excess chemical potential. 
There is overall quite good agreement, although the weakening of the RHNC 
results with increasing density and decreasing temperature, a shortcoming 
common to all integral equation closures, does become evident. 

In tables 1 and 2 we report a few significant comparisons between RHNC 
thermodynamic and structural results and recent numerical simulation data 
[1^ [H] for SW systems of different range. 

We notice a progressive worsening of thermodynamic results for stronger 
couplings, meaning lower temperatures and higher densities, and for decreas- 
ing A. Such behavior as a function of well width is understandable, since in 




(21) 
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/5P/p 




0.5 1 1.5 2 



/5e 

Figure 3: SP-MC pressure vs. inverse temperature for a SW fluid with 
A = 1.5 obtained by Labfk et al. [4Dj at p* = 0.8 (filled circles) and p* = 0.9 
(empty circles). Solid lines are the RHNC results. 
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Figure 4: SP-MC chemical potential vs. inverse temperature for a SW fluid 
with A = 1.5 obtained by Labik et al. [10] at p* = 0.8 (flUed circles) and 
p* = 0.9 (empty circles). Solid lines are the RHNC results. 
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Table 1: Comparison of simulation and RHNC thermodynamic quantities of 
square-well fluids. For each density, the first row contains Monte Carlo data 
of Largo and Solana [12] and the second row our results. 



A = 1.2 p* 
0.1 

0.8 

A = 1.5 

0.1 
0.7 

A = 2.0 

0.1 
0.6 
0.7 



;P/p -U/Ne 

T* = 0.7 

0.76 0.62 
0.76 0.62 

0.77 3.68 
0.68 3.58 

T* = 1.5 

0.78 0.92 
0.78 0.92 

2.31 5.28 
2.30 5.27 

T* = 3.0 

0.69 1.96 
0.69 1.98 

1.19 9.75 
1.18 9.76 

2.01 11.35 
1.97 11.35 



(3P/p -U/Ne 

T* = 1.0 

0.96 0.40 
0.96 0.40 

2.72 3.43 
2.68 3.38 

T* = 3.0 

1.03 0.69 
1.03 0.69 

4.02 5.08 
3.98 5.07 

T* = 5.0 

0.91 1.72 
0.92 1.73 

2.44 9.58 
2.43 9.59 

3.51 11.16 
3.47 11.16 



(3P/p -U/Ne 

T* = 3.0 

1.17 0.23 
1.17 0.23 

6.03 3.09 
5.99 3.09 
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Table 2: Comparison of simulation and RHNC contact values g{o'~^), g{Xa~) 
and g{X(7~^) of the radial distribution function of square-well fluids. Monte 
Carlo data (MC) from Largo et al. [8j. 



g{a+) giXa-) g{\a+) 



X T* p* MC RHNC MC RHNC MC RHNC 



1.1 


0, 


.5 


0, 


.1 


7.254 


7.260 


7.179 


7.179 


0, 


.966 


0, 


.973 








0, 


.5 


6.088 


5.835 


5, 


,816 


5.460 


0, 


.787 


0, 


,739 








0, 


.7 


5.667 


5.681 


5.307 


4.967 


0, 


.718 


0, 


,672 


1.2 


0, 


.7 


0, 


.1 


4.126 


4.150 


4.030 


4.046 


0, 


.968 


0, 


.970 











.5 


3.449 


3.272 


3, 


,122 


2 


.945 


0, 


.748 


0, 


,706 








0, 


.7 


3.281 


3.186 


2, 


,810 


2 


.744 


0, 


.675 


0, 


,658 


1.5 


1, 


.5 





.1 


1.952 


1.957 


1, 


,832 


1 


.835 


0, 


.941 


0, 


.944 








0, 


.5 


1.989 


1.994 


1, 


,382 


1 


.382 


0, 


.709 


0, 


,709 











.7 


2.783 


2.770 


1, 


,115 


1 


.149 





.590 


0, 


,590 




2, 


.0 


0, 


.1 


1.661 


1.655 


1, 


,563 


1 


.561 


0, 


.945 


0, 


.947 








0, 


.5 


2.006 


2.003 


1, 


,274 


1 


.272 


0, 


.771 


0, 


,772 








0, 


.7 


2.880 


2.865 


1, 


,063 


1 


.061 


0, 


.646 


0, 


,644 


1.8 


2, 


.0 


0, 


.1 


1.839 


1.914 


1, 


,599 


1 


.642 


0, 


.974 


0, 


.996 








0, 


.5 


2.199 


2.211 


1, 


,175 


1 


.181 


0, 


.714 


0, 


,716 








0, 


.7 


3.487 


3.474 


1. 


,141 


1 


.147 


0, 


.693 


0, 


,695 




3, 


.0 


0, 


.1 


1.449 


1.456 


1, 


,320 


1 


.320 


0, 


.942 


0, 


.945 








0, 


.5 


2.162 


2.162 


1, 


,079 


1 


.080 


0, 


.773 


0, 


,774 








0, 


.7 


3.392 


3.365 


1. 


,043 


1, 


.045 


0, 


.747 


0. 


,751 
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the limit of very large A and vanishing e, the properties of SW fluid should 
be well described by HS physics, supplemented by the van der Waals mean 
field approximation. In such a limit RHNC should provide very accurate 
results. In the opposite case, we would expect the maximum deviation from 
the hypothesis of universality of the bridge function and then the maximum 
discrepancy from simulation. In the only previous RHNC study of SW fluid 
that we are aware of, ref. [I6], a limited analysis of the differences between 
HS and SW bridge functions has been attempted. In that paper, the authors 
assumed the same functional form for the bridge function outside the hard 
core for both HS and SW systems and determined the parameters via a mean- 
square fitting of Monte Carlo (MC) data via RHNC. Apparently, optimal SW 
bridge functions should be slightly smaller than the best HS functions for the 
same system, but we note that the constrained procedure did not allow for 
qualitative changes of the functional form. We can add that the most im- 
portant for possible changes is the region within and close to the hard core. 
Direct evidence for this comes from a few calculations we performed by mod- 
ifying the long range part of the bridge function in a very crude way: we 
made it vanish starting from its first zero. The resulting RHNC solution was 
almost unaffected by this change, thus confirming that any future effort to 
improve the bridge function beyond the HS approximation should be focused 
on the core interval between zero and the HS diameter. 

We have studied the RHNC liquid- vapor coexistence of four SW systems 
of width size A = 1.25, 1.5, 1.75, and 2. In figures 5-8 we present our 
results for the coexistence curves of the four SW fluids we have examined. 
RHNC data (squares connected by a line) are compared with different sets 
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Figure 5: Vapor-liquid T* versus p* coexistence for the SW fluid of range 
A = 1.25: squares, this work (circled squares indicate extrapolation of the 
vapor branch), a line has been drawn through the points as a visual guide; 
dots, Vega et al. p]; crosses, del Rio et al. [5]. 
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Figure 6: Vapor-liquid T* versus p* coexistence for the SW fluid of range 
A = 1.5: squares, this work (circled squares indicate extrapolation of the 
vapor branch), a line has been drawn through the points as a visualguide; 
dots, Vega et al. p]; crosses, del Rio et al. [5]. 
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Figure 7: Vapor-liquid T* versus p* coexistence for the SW fluid of range 
A = 1.75: squares, this work (circled squares indicate extrapolation of the 
vapor branch), a line has been drawn through the points as a visualguide); 
dots, Vega et al. p]; crosses, del Rio et al. [5]. 
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Figure 8: Vapor-liquid T* versus p* coexistence for the SW fluid of range 
A = 2.0: squares, this work (circled squares indicate extrapolation of the 
vapor branch), a line has been drawn through the points as a visualguide; 
dots, Vega et al. p]; crosses, del Rio et al. [5]; triangles, de Miguel [4j. 
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of simulation data (see figure captions). All the cases where RHNC did not 
provide a real crossing between the two branches (low and high density) of 
the f3P-(3fi curves have been marked by a circle around the square. 

Two main features of the results are evident: the overall semiquantita- 
tive accuracy of the coexistence curve in all the investigated cases and the 
progressive worsening of the quality of the RHNC results as the range A 
is reduced. The latter appears as the obvious consequence, at the level of a 
phase diagram, of our observations about the quality of thermodynamic data 
as a function of the range. At A = 2, our results closely follow the data of de 
Miguel [3] and of del Ri'o et al. [5j , while Gibbs ensemble Monte Carlo results 
by Vega at al. [3] seem to be definitely biased toward higher densities. In 
addition, the unphysical changes of the curvature seen in the MC data sug- 
gest that one might consider quite optimistic the published statistical error 
bars. 

At A = 1.75, RHNC results start to show clear evidence of a distortion 
of the coexistence curve at the highest temperature, although all the points 
below T* = 1.70 are in good agreement with MC data. Since the extrapola- 
tion at this temperature corresponds to the case shown in figure 2, while the 
remaining extrapolations correspond to closer missed intersections, we have 
a measure of the quantitative effect of our extrapolation procedure on the 
resulting binodal line. 

Data at A = 1.5 and A = 1.25 clearly show a progressive decrease of the 
quality on the high density side of the curve. As previously discussed, this 
region corresponds to the region where RHNC could be substantially affected 
by improvements in the description of the bridge function. We stress that 
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coexistence is a severe test for the quality of integral equations. In particular, 
pressure and chemical potential are quantities quite sensitive to small changes 
of the closure. 

However, according to Liu et al. [9j, at A = 1.25 the liquid- vapor transi- 
tion becomes metastable and is preempted by the freezing transition. Thus, 
the observed decreased quality at smaller widths is somewhat compensated 
by the fact that the integral equation results for such case are referring to a 
metastable liquid. 

In all the cases, integral equations do not allow one to get close enough 
to the critical point to provide a direct estimate of its location. However, 
a good quality in reproducing the low temperature part of the liquid vapor 
coexistence and the hypothesis of Ising-like value for the critical exponent 
(/? = 0.325) may provide a first approximation for the extent of the two- 
phase region. 

Conclusions 

We have shown that RHNC is able to reproduce the low-temperature part of 
the liquid-vapor phase diagrams of square-well fluids. With increasing tem- 
perature and depending on the extent of the attractive well, there is a region 
of the binodal line that can be obtained by a slight extrapolation of the vapor 
phase beyond the numerical limit of the quasi- spinodal line. Eventually, ap- 
proaching the critical temperature, the extrapolation becomes very poor and 
the shape of the coexistence line is strongly affected. The agreement with 
the best computer simulations is good for the widest well sizes, but becomes 
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less satisfactory for A at or below 1.5. Since the only approximation of the 
theory is the bridge function (apart from inherent limitations of numerical 
computation), our results point to the need for an improvement of the short 
range part of the bridge function. 

It is interesting to note also that our results are consistent with similar 
RHNC calculations on the Lennard- Jones fluid fTll [2S] , where it is possible 
to obtain a larger part of the binodal line than in the case of A = 2, thus 
confirming our findings on better performances for wider range potentials. 

We stress again that the main purpose of the present investigation is 
not to propose RHNC as an alternative to more refined theories such as 
SCOZA or HRT when a highly accurate description of the critical properties 
is required. Rather, we have the more limited goal of checking the ability of 
the RHNC formulation to provide an approximate but quick tool for locating 
the liquid-vapor phase transition in SW-like systems. In all the isotropic 
SW cases we have studied, the theory allows one to predict from scratch 
an unknown phase coexistence in a computational time that is negligible 
compared to numerical simulations, using existing numerical techniques for 
integral equation solutions. Such a feature may only be marginally interesting 
for isotropic SW fluids, but becomes crucial for anisotropic patchy models of 
globular macromolecules. We are currently investigating this topic. 

After this paper was submitted for publication, a very thorough work 
on the same topic of liquid-vapor phase equilibrium for the SW fluid us- 
ing a different integral equation was published by El Mendoub et al. [13] . 
These authors use an efficient "adaptive grid" 01] that automates the map- 
ping of the no- solution space in the temperature-density plane. Comparing 
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their results to ours, we notice that RHNC results for the radial distribution 
function show a uniform better agreement with simulation data. Since de- 
viations from computer simulation binodal lines shown by their closure and 
the present study have opposite directions, a detailed comparison between 
their closure and RHNC would be interesting . 
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